Package flip

su CRAN ma soprattutto su github

To install this github version type (in R):

library(devtools)

install_github("flip", "livioivil")

(if devtools is not installed yet: install.packages("devtools") )

seeds data

Leggo i dati seeds

library(flip)
data(seeds)

seeds
##    grp    x     y
## 1    0   NA    NA
## 2    0   NA    NA
## 3    0   NA    NA
## 4    0   NA    NA
## 5    0   NA    NA
## 6    0   NA    NA
## 7    0   NA    NA
## 8    0   NA    NA
## 9    0 6.03 12.54
## 10   0 4.20 14.81
## 11   0 4.49 16.71
## 12   0 2.00  7.53
## 13   0 2.84  7.02
## 14   0 3.88  8.09
## 15   0 2.04  5.76
## 16   0 5.48 18.01
## 17   0 2.31  8.81
## 18   0 1.90  8.17
## 19   0 1.75  6.62
## 20   0 3.02  7.69
## 21   1   NA    NA
## 22   1   NA    NA
## 23   1   NA    NA
## 24   1 3.31 18.49
## 25   1 6.56 19.20
## 26   1 3.16  9.85
## 27   1 4.07 15.83
## 28   1 2.09  6.16
## 29   1 6.72 17.58
## 30   1 3.93 19.29
## 31   1 2.56 10.77
## 32   1 8.30 18.31
## 33   1 4.21 10.56
## 34   1 1.86  9.48
## 35   1 3.09 12.54
## 36   1 5.09 18.35
## 37   1 4.08 11.84
## 38   1 3.63 11.44
## 39   1 2.61  7.66
## 40   1 5.21 12.00
# leave out the NAs
dati=seeds[!is.na(seeds$x),]

Analisi descrittive:

#funzione che calcola qualche statistica descrittiva
stats=function(v){
  sds=apply(v,2,sd)
  sds=round(sds,3)
  sds=paste("sd     :",sds,sep="")
  rbind(summary(v),sds)}

#la applico separatamente nei due campioni
by(dati[,2:3],dati$grp,stats)
## dati$grp: 0
##           x                 y           
##     "Min.   :1.750  " "Min.   : 5.760  "
##     "1st Qu.:2.030  " "1st Qu.: 7.402  "
##     "Median :2.930  " "Median : 8.130  "
##     "Mean   :3.328  " "Mean   :10.147  "
##     "3rd Qu.:4.272  " "3rd Qu.:13.107  "
##     "Max.   :6.030  " "Max.   :18.010  "
## sds "sd     :1.467"   "sd     :4.228"   
## -------------------------------------------------------- 
## dati$grp: 1
##           x                 y          
##     "Min.   :1.860  " "Min.   : 6.16  "
##     "1st Qu.:3.090  " "1st Qu.:10.56  "
##     "Median :3.930  " "Median :12.00  "
##     "Mean   :4.146  " "Mean   :13.49  "
##     "3rd Qu.:5.090  " "3rd Qu.:18.31  "
##     "Max.   :8.300  " "Max.   :19.29  "
## sds "sd     :1.753"   "sd     :4.354"

box-plot

par(mfrow=c(1,2))
boxplot(seeds$x~seeds$grp,main="x",col=c(0,2))
boxplot(seeds$y~seeds$grp,main="y",col=c(0,2))

Verifica di ipotesi

Confronto:

per entrambe le variabili: z=x,y

res=flip(x+y~grp,data=dati,perms=5000)
##     
summary(res)
##  Call:
##  flip(Y = x + y ~ grp, data = dati, perms = 5000) 
## 4999 permutations.
##   Test  Stat tail p-value sig.
## x    t 1.320   ><  0.1920     
## y    t 2.061   ><  0.0532

In questo caso è ragionevole fare ipotesi alternative unidirezionali

H1(z): F(z|grp=0)<F(z|grp=1)

res=flip(x+y~grp,data=dati,perms=5000,tail=1)
##     
summary(res)
##  Call:
##  flip(Y = x + y ~ grp, data = dati, tail = 1, perms = 5000) 
## 4999 permutations.
##   Test  Stat tail p-value sig.
## x    t 1.320    >  0.0950     
## y    t 2.061    >  0.0202    *

Visualizziamo gli istogrammi delle distribuzioni delle statistiche test di permutazione (calcolate sull’orbita dei dati osservati).

##set the following if you get an error (mostly using Rstudio)
#par(mar=c(1,1,1,1))
hist(res)

E’ però fondamentale studiare la distribuzione congiunta delle due statistiche test:

plot(dati$x,dati$y,pch=21,bg=dati$grp+1)
#medie di gruppo:
XYbars=rbind(colMeans(dati[dati$grp==0,2:3]),colMeans(dati[dati$grp==1,2:3]))

points(XYbars[,1],XYbars[,2],pch=21,bg=c(1:2),cex=2)
abline(h=mean(dati$y))
abline(v=mean(dati$x))

La dipendenza tra le due variabili induce una dipendenza anche nella distribuzione congiunta delle statistiche test

plot(res) 
abline(h=0)
abline(v=0)

E di conseguenza anche nei corrispettivi valori lambda (i p-value)

#calcolo i lambda per ogni permutazione, per ogni variabile
res@permP=t2p(res)
#scatter plot:
plot(res@permP[,1],res@permP[,2],col="#F98400",bg="#F2AD00",pch=21)
points(res@permP[1,1],res@permP[1,2],bg="#00A08A",col="#F2AD00",cex=2,pch=21)
abline(v=p.value(res)[1])
abline(h=p.value(res)[2])

ESERCIZIO

è opportuno ricordare che la dipendenza tra statistiche test nella distribuzione congiunta non altera le distribuzioni marginali (cioè univariate) delle statistiche test.

Farte una verifica empirica con comandi del tipo hist(res@permP[,"x"]) o meglio plot(ecdf(res@permP[,"x"]))

Combinazione Non Parametrica di test

Confronto le ipotesi:

Funzione di Fisher

(res.fish=npc(res,comb.funct = "F")) #sono implementate molte altre funzioni di combinazione e opzioni, si veda ?npc
##    comb.funct nVar  Stat p-value
## V1     Fisher    2 6.256  0.0374

Visualizziamo

# calcolo il vettore dei p-values sulla base del vettore di statistiche test
res.fish@permP=t2p(res.fish)

#vettore logico che individua tutti i p-value (ottenuti da permutazioni casuali) che sono <= al p-value combinato osservato
lower=res.fish@permP<=res.fish@permP[1]

#plot
plot(res@permP[,1],res@permP[,2],col="#F2AD00",bg=c("#F98400","#FF0000")[lower+1],pch=21,main="p-values",asp=1)
points(res@permP[1,1],res@permP[1,2],bg="#00A08A",col="#F2AD00",cex=2,pch=21)

#la statistica test osservata è pari a -6.19138 = -log(p-value x) -log(p-value y) 
# ( salvato in res.fish@permP[1] )
#la seguente iperbole è il luogo dei punti con la stessa statistica test:
curve(exp(-6.19138)/x,ylim=c(0,1),from=0.000001,to=1,add=TRUE)

Visualizziamo i corrispettivi punti nello spazio delle statistiche test

plot(res,bg=c("#F98400","#FF0000")[(lower)+1],asp=1)

Funzione di Tippett (min-p)

(res.tip=npc(res,comb.funct = "minP")) 
##    comb.funct nVar Stat p-value
## V1       minP    2 49.5  0.0342

Visualizziamo

res.tip@permP=t2p(res.tip)
lower=res.tip@permP<=res.tip@permP[1]

plot(res@permP[,1],res@permP[,2],col="#F2AD00",bg=c("#F98400","#FF0000")[lower+1],pch=21,main="p-values",asp=1)
points(res@permP[1,1],res@permP[1,2],bg="#00A08A",col="#F2AD00",cex=2,pch=21)
abline(h=min(p.value(res)))
abline(v=min(p.value(res)))

Controllo della molteplicità

Errore di I tipo per variabile x : P(p(x)<=alpha|H0)

Nei test esatti questa probabilità è <=alpha

Altrettanto vale per y.

FamilyWise Error Rate (FWER) - Errore di I tipo multiplo

Qual è a probabilità di commettere ALMENO un errore facendo i due test?

P(p(x)<=alpha OR p(y)<=alpha|H0)=?

NON conosciamo la vera configurazione di H0

Non sappiamo cioè se:

  1. H0(x) AND H0(y) : si possono fare 0, 1 o 2 errori
  2. H0(x) AND H1(y) : si possono fare 0 o 1 errore
  3. H1(x) AND H0(y) : si possono fare 0 o 1 errore
  4. H1(x) AND H1(y) : nessun errore possibile

e quindi 1. H0(x) AND H0(y) : test combinato 2. H0(x) AND H1(y) : test su x (univariato) 3. H1(x) AND H0(y) : test su y (univariato) 4. H1(x) AND H1(y) : nessun test

Bonferroni

Combinazione di Fisher

plot(res@permP[,1],res@permP[,2],col="#F2AD00",bg="#F98400",pch=21,asp=1)
points(res@permP[1,1],res@permP[1,2],bg="#00A08A",col="#F2AD00",cex=2,pch=21)

# p-value combinato : (caso I)
curve(exp(-6.19138)/x,ylim=c(0,1),from=0.00001,to=1,add=TRUE,col="red",lwd=2)

# p-value per H0(x): (caso II)
abline(v=p.value(res)[1],col="green",lwd=2)

# p-value per H0(y): (caso III)
abline(h=p.value(res)[2],col="blue",lwd=2)

legend("topright",legend=c("H0(x,y) (Fisher)","H0(x)","H0(y)"),col=2:4,bty="n",lwd=2)

Comb. di Tippett (min-p)

plot(res@permP[,1],res@permP[,2],col="#F2AD00",bg="#F98400",pch=21,asp=1)
points(res@permP[1,1],res@permP[1,2],bg="#00A08A",col="#F2AD00",cex=2,pch=21)

# p-value combinato : (caso I)
abline(v=min(p.value(res)),col="red",lwd=2)
abline(h=min(p.value(res)),col="red",lwd=2)

# p-value per H0(x): (caso II)
abline(v=p.value(res)[1],col="green",lwd=2)

# p-value per H0(y): (caso III)
abline(h=p.value(res)[2],col="blue",lwd=2,lty=2)

legend("topright",legend=c("H0(x,y) (Tippett)","H0(x)","H0(y)"),col=2:4,bty="n",lwd=2)

Correzione dei p-value via Fisher (lento quando ci sono molte variabili)

(res=flip.adjust(res,"F"))
##   Test  Stat tail p-value Adjust:Fisher
## x    t 1.320    >  0.0950        0.0950
## y    t 2.061    >  0.0202        0.0374

Correzione dei p-value via minP

(res=flip.adjust(res,"minP"))
##   Test  Stat tail p-value Adjust:Fisher Adjust:minP
## x    t 1.320    >  0.0950        0.0950      0.0950
## y    t 2.061    >  0.0202        0.0374      0.0342

Regioni di rifiuto

Combinazione di Fisher

Combinazione di Tippett (min-p)